#ifndef EDS_HAMILTONIAN_ACTOR_HEADER
#define EDS_HAMILTONIAN_ACTOR_HEADER

#include<pthread.h>
#include<vector>

#include<src/matrix_actor.h>
#include<src/configuration.h>

namespace eds {

  template <class SCALAR> void * thread_operate_by_A(void * arg);

  template <class SCALAR> class HamiltonianActor : public pcps::MatrixActor<SCALAR> {

    private:

      /// \brief  the number of threads to use
      int _nthread;

      /// \brief  the number of sites
      int _ns;

      /// \brief  the total occupation of the lattice
      int _to;

      /// \brief  the number of configurations on each site
      int _cps;

      /// \brief  threshold for convergence
      const double _thresh;

      /// \brief  offests used in index computation
      std::vector<int> _offsets;

      /// \brief  spacings used in index computation
      std::vector<int> _spacing;

      /// \brief  z array lookup table used in index computation
      std::vector<int> _z_lookup;

      /// \brief  z arrays table used in index computation
      std::vector<int> _z_array;

      /// \brief  integer workspace used in index computation
      std::vector<int> _work;

      /// \brief  coefficients associated with spin ladder raising operator
      std::vector<double> _plus_factor;

      /// \brief  coefficients associated with spin ladder lowering operator
      std::vector<double> _minus_factor;

      /// \brief  bond indices
      std::vector<int> _bonds;

      /// \brief  bond operator coefficients
      std::vector<double> _bond_strengths;

      /// \brief  triangle indices
      std::vector<int> _triangles;

      /// \brief  lattice configuration
      std::vector<eds::Configuration> _start_configs;

      /// \brief  compound indices of each thread's end configuration
      std::vector<int> _end_cmpds;

      /// \brief  space for the threads to store their multiplication results
      std::vector<SCALAR> _thread_workspace;

      /// \brief  pointers to each thread's storage space
      std::vector<SCALAR *> _thread_work_ptrs;

    public:

      // constructor
      HamiltonianActor(const int ns,
                       const int to,
                       const int cps,
                       const std::vector<int> & bonds,
                       const std::vector<double> & bond_strengths,
                       const std::vector<int> & triangles,
                       const int nthread);

      // function to perform y = A d
      void operate_by_A(const int n, SCALAR * const d, SCALAR * const y, const int myrank, const int tid);

      // function to perform y = M d
      void operate_by_M(const int n, SCALAR * const d, SCALAR * const y, const int myrank, const int tid);

      // function to perform y = M^(-1) d
      void operate_by_M_inv(const int n, SCALAR * const d, SCALAR * const y, const int myrank, const int tid);

      // function to check if the iterations have converged
      // result should be computed using the root process/thread and returned on all processes/threads
      bool converged(const double residual, const int myrank, const int tid);

      void apply_bond(eds::Configuration & config,
                      const int p,
                      const int q,
                      const int curr_cmpd,
                      const SCALAR c,
                      SCALAR * const y,
                      int * const work);

      //void triple_products(const int n, const SCALAR * const x);
      void triple_products(const int n, const int ntpv, const int tri_num);

      SCALAR tp_term(const int i, const int j, const int k, const SCALAR * const x,
                     eds::Configuration & config);

      struct ThreadArgs {
        int tid;
        int myrank;
        eds::HamiltonianActor<SCALAR> * obj_ptr;
        pthread_mutex_t * mtx_ptr;
        int n;
        const SCALAR * d;
        SCALAR * y;
      };

      friend void * eds::thread_operate_by_A<SCALAR>(void *);

  };

}

#endif
